

PMaxEst = function(ratio, gop, K){
    
    ### Step 2: get ratio against the "baseline" p_{0...0}
    d1 = 2^(2*K + 1)
    ratio.base = numeric(d1)
    ratio.base[1] = 1
    ratio.base[2] = ratio[1,1] * ratio.base[1]
    ratio.base[3:4] = ratio[2,1:2] * ratio.base[1:2]
    ratio.base[5:8] = ratio[3,] * ratio.base[1:4]
    if(K > 2) for(k in 2:(K-1)){
        # change lk
        ratio.base[(2^(2*k-1)+1) : (2^(2*k))] = rep(ratio["lk",],
                                                    each = 2^(2*k-3)) * 
            ratio.base[1: (2^(2*k-1))]
        # change ak
        ratio.base[(2^(2*k)+1) : (2^(2*k+1))] = rep(ratio["ak",],
                                                    each = 2^(2*k-2)) * 
            ratio.base[1: (2^(2*k))]
    }
    ratio.base[(2^(2*K-1)+1) : (2^(2*K))] = rep(ratio["lK",],
                                                each = 2^(2*K-3)) * 
        ratio.base[1: (2^(2*K-1))]
    ratio.base[(2^(2*K)+1) : (2^(2*K+1))] = rep(ratio["aK",],
                                                each = 2^(2*K-2)) * 
        ratio.base[1: (2^(2*K))]
    k = ratio.base / max(ratio.base)
    
    f = function(x){
        # prod(k) * (x^d1) / prod(1-x*k) - gop
        sum(log(k)) + d1 * log(x) - sum(log(1-x*k)) - log(gop)
    }
    
    tol = 10
    if(f(1-0.1^tol)<0) p.max = 1 else{
        p.max = uniroot(f, lower=0, upper=1, tol=10^(-tol))$root    
    }
    
    p = k * p.max
    
    return(list(p.max = p.max, p = p))
    
}




